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Abstract 

A large class of physically important nonlinear and nonhomogeneous evolution problems, charac- 
terized by advection-like and diffusion-like processes, can be usefully studied by a time-differential 
form of Kolmogorov's solution of the backward-time Fokker-Planck equation. The differential so- 
lution embodies an integral representation theorem by which any physical or mathematical entity 
satisfying a generalized nonhomogeneous advection-diffusion equation can be calculated incremen- 
tally in time. The utility of the approach for tackling nonlinear problems is illustrated via solution 
of the noise-free Burgers and related Kardar-Parisi-Zhang (KPZ) equations where it is shown that 
the differential Kolmogorov solution encompasses, and allows derivation of, the classical Cole-Hopf 
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and KPZ transformations and solutions. A second example, illustrating application of this ap- 
proach to nonhomogeneous evolution problems, derives the Feynman-Kac formula appropriate to 
a Schrodinger-like equation. 

Key words, nonlinear partial differential equations, stochastic methods, Cole-Hopf transforma- 
tion, Feynman-Kac formula 
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1. Introduction 

A wide array of physically important time-evolutionary problems, in areas as diverse as hydrody- 
namic turbulence [1-3], quantum mechanics [4,5], and the evolution of cosmic large-scale struc- 
ture [6], can be modeled by nonlinear equations of the generic form 

ije + b • Vi] + zA7 2 r/ = f(r), //,.//,,...„; (1) 

where rj represents a scalar field whose evolution within a finite or non-finite space-time domain is 
determined, at least in part, by processes that in some fashion, are analogous to advection and/or 
molecular diffusion within a hydrodynamic flow. Here, b and v correspond to, or are analogs 
of, the local dimensional or nondimensional drift field and diffusion coefficient, respectively, while 
the source term, /, can include linear and nonlinear terms in both r\ and derivatives of r] up to 
order n. In addition, since we will conceptualize the incremental, stochastic construction of r\ as 
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a procedure in which random walkers are launched from a space-time solution point (x, 6) toward 
known/previously calculated r/, we will follow Kolmogorov and state the evolution of r\ in terms of 
backward time, 6. 

It is well known, particularly in the applied mathematics community [7, 8] that linear models, 
which are encompassed by ([I]), can be solved over finite (non- incremental) time intervals via appli- 
cation of Kolmogorov's solution of the backward-time Fokker Planck (FP) equation. In cases where 
(PQ) has either a nonlinear or linear source term, /, in rj, or a drift b and/or diffusivity dependent 
on rj (where the latter can arise in more general forms of the third term on the left of ([I])), a 
Kolmogorov solution can still be expressed in representative or schematic form [9,10]; however, in 
these instances, one or more stochastic path integrals in the a priori unknown scalar field, rj(x,6), 
appear. While one can, in principal, formulate a system of equations in the unknown rj, and attempt 
a numerical solution over the entire space-time domain, Q, this is computationally expensive and 
apparently remains untried. Here, we focus on simpler, low cost, and often analytically tractable 
local solutions. 

The purpose of this paper is to point up the utility of applying a time- differential version 
of Kolmogorov's solution to nonlinear and nonhomogeneous evolution problems of the form in 
(PQ). In order to demonstrate the potential of this approach for tackling nonlinear problems, we 
show how the differential Kolmogorov solution can be used to derive the Cole-Hopf and related 
Kardar-Parisi-Zhang (KPZ) solutions of the noise-free Burger's [11, 12] and KPZ [13] equations. 
We then apply this approach to nonhomogeneous problems, deriving the Feynman-Kac formula as 
a solution to a Schrodinger-like equation. We find that, at least in these test cases, the differential 
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Kolmogorov solution offers an alternative, relatively straightforward, physically transparent basis 
for solution. Moreover, since this approach can be adapted to accommodate the effects of boundary 
interactions/boundary conditions on the evolution of r/ [14] (see also the caption to Fig. 2), the 
potential for further applications seems clear. 

Our jumping-off point is thus Kolmogorov's solution of the backward-in-time Fokker-Planck 
equation, connecting the stochastic kinematics of random walkers within a finite or non-finite 
space-time region, Q = D x (9,T], to the evolution of the scalar field, r/(x,9), within Q. Here, D 
is the spatial domain, 9 = T — t, is again backward time, t forward time, and T the forward time 
at which a solution for ry(x, 9) is sought. The solution, written here for finite Q, and represented 
schematically in Fig. 1, reads [7] : 

77(x,0) =E X)S [«7(x(r),r)] + E ^{ X {r))] - E g[ f f{ X {s),s)ds] (2) 

Je 

where the first two terms on the right capture the effects of the random walker (RW) swarm, 
launched from the solution point (x, 9), sampling (typically known) Dirichlet conditions, r/ = 
g(x, #'), on the boundary 5Q of Q, and the final condition r] = (j)(x,9 = T) on the final time- 
slice, D x {8 = T}. In cases where Neumann and/or mixed conditions are imposed on SQ, an 
additional term similar to the first appears [14]. Throughout, E will represent the expectation 
taken with respect to the solution point. In addition, in ([2]), r represents the random time at which 
any given RW impinges either on 5Q (first term on right side of ([2|)), or on D x {9 = T} (second 
term) . 

Couching the development mainly in hydrodynamic terms, and depicting the time-incremental 
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Figure 1: Kolmogorov's local, non-time-incremental stochastic construction of the scalar ?y(x, #), 
as given by Eq. ([2]), appropriate to linear versions of Eq. ([1]), where the drift, b, and diffusivity, v 
(or their non-hydrodynamic analogs) are known. In this illustration, r) is subject to the Dirichlet 
condition r] = g(x, 9) on the boundary 5Q of the space-time cylinder Q, while on the final-time 
slice, D x {9 = T}, r/ = 0(x). 

solution in Fig. 2, we focus below on a differential form of ([2]), expressed over short backward time 
intervals, A9 = Oj+i — 9j, where A9 is chosen so that the magnitudes of all three terms on the 
left side of (pQ) remain comparable over the interval. This is an important point both physically 
and computationally and translates to the condition |(x — 5L{9jj r \))\ 2 /l 2 D ~ [b(x, 9j)A9] 2 /(uA9) = 
blb/u = O(l), where b = |b(x, 9)\, h(x,9j) = —u(x,9j) is the local drift, Id = V vA9 is the local 
diffusion length scale, lb = bA9 is the local advection length scale, and x(#j + i) is the mean position 

* 

of the random walk swarm at Oj+i- [Note, 5c(#j+i) = x + E h(x(s')ds' = x + b(x, 9j)A6 + 



Thus, in hydrodynamic terms, we arrive at the requirement that A9 must be chosen so that 
the local Reynolds number, Re = bl^/v, remains order 1, or more meaningfully, so that the local 



0(A9 2 )}. 
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Figure 2: Incremental stochastic construction of rj(x.,9j), appropriate to nonlinear versions of Eq. 
([T]). Here, an approximate 'cone-of-influence' is shown, indicating how the local solution at (x, Oj) 
would typically depend on previously determined 77, computed incrementally over [Oj+i,T). In this 
illustration, a drift field (or drift field analog), b, having a significant x\ component, pushes RW's 
predominantly to the right. Near Dirichlet boundaries, RW impingement with a boundary element 
$Qk = SD x A9k (where 5D is the boundary of D and A9k = 9k+i — 9k), is readily accommodated 
using the first expectation in Eq. ([2]), applied over 5Qk- A similar (or additional) term would be 
required near Neumann (or mixed) boundaries [14]. 

advection length scale, If, = bA9, is on the order of the local diffusion length scale, Id- More 
importantly, this condition implies that non-fixed time steps may, under various circumstances, be 
required in both analytical [14] and Monte Carlo solutions, a point recognized in the applied math 
literature [15], but apparently not well known in the engineering and physics communities. 

In order to allow comparisons with known results, we will: i) assume that the solution point 
(x, 9) £ Q is sufficiently removed from boundaries to allow neglect of boundary interaction effects, 
and ii) focus on one-dimensional problems. The method to be described carries over without 
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difficulty to multiple dimensions, and as noted, can be adapted to handle boundary conditions. In 
addition, for notational simplicity, we will set v = 1; based on the discussion immediately above, 
this corresponds to nondimensionalizing ([T]) using If, (or equivalently, Id), b, and A9 for the length, 
drift, and time scales, respectively, and then to setting Ret, = 1. 

Thus, we express the time-differential form of © in either of two equivalent forms: 



V(x,9 j )= I p{x,dj\x\Bj +1 ))ri{x' ,6 j+1 )dx' - E \ ^ f( X (s'),s')ds' (3) 

-oo 



riWi) - E [ n{Xi ^ +l) ]v(4,0 j+ i) - ^ f(Xk(s%s')ds' (4) 



or 

fe=l 

where p is the transition density, n(x[, 9j+i)/N is the fraction of N RW's launched from the solution 
point (x,9j) reaching a line segment, Axi (on the x-axis), during the interval (0j,9j+i\, x i £ Axi, 
and N x is the number of increments used to discretize a length 2L of the x-axis. Here, 2L is, e.g., 
centered on x, and L is chosen to be much larger than the maximum RW displacement produced 
by combined advection and diffusion over A9. 

With regard to the differential solutions in ([3|) and we note that ([3]) is useful when the 
transition density, p, can be determined analytically via the FP equation; this approach is illus- 
trated in our derivation of the Cole-Hopf/KPZ (CH/KPZ) transformation below. The second 
form, (U]), which is useful, for example, in Monte Carlo simulations where analytical p's are not 
available, is based on the following simple relationship between p and the fraction n^x^Oj+^/N : 
p(x,9j\x i ,9j+i) = limAzi-^o n (^i) ^i+i)/ (NAxi). As a practical aside, and considering the frac- 
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tions n(x i ,9j-^-i)/N are obtained by launching a swarm of N RW's from the solution point (x,9j), 
and then by counting the number of RW's that subsequently reach each Axi over the increment 
(0j,9j+i\. 

With regard to the path integral terms in ([3]) and ([3D, these can be expressed, depending on 
the application, as f(x, 6j + i)A8 + 0(A9 2 ), or f(x, 9j)A0 + 0(A8 2 ). This point is exploited in our 
derivation of the Feynman-Kac formula below. 

From a mathematical standpoint, the power of both the incremental and non-incremental Kol- 
mogorov solution, ([3]) (or (J3|)) and ((2|), respectively, derives from the fact that these embody repre- 
sentation theorems [see, e.g., [9]] by which, in the present context, any entity, /J(p, a>), satisfying a 
generalized transport equation of the form f} u + v • V p /3 + = F((3), and having an associated 

stochastic differential equation (SDE), dx(uS) = v(p, u) + \^2Kdw(oj), can be determined in time- 
incremental fashion. Here, p and u represent variables that are, in the sense of this generalized 
transport equation, space- and time-like, respectively, v and k are, respectively, drift field and 
diffusion coefficient analogs, F is a generalized source term, and w is a Wiener process. 

From a physical standpoint, and (J3|) likewise embody powerful features. In particular, since 
solutions of ([I]) via ([3]), or equivalently ([3]), require determination of the transition density p (or 
equivalently, n(x[, 9j + \) / (N Ax^)) , where p typically captures short time- and length-scale physics, 
then in constructing solutions, we are immediately forced to confront physical features extant 
on these scales. In many continuum-scale scalar transport problems, for example, one can often 
reasonably assume that small scale statistical features are adequately captured by the Fokker- 
Planck, or somewhat more generally, the differential Chapman-Kolmogorov (CK) equation [16]. In 
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these cases, and again as illustrated in the first example below, due to the local nature of ([3]) or 
©) P typically corresponds to a fundamental solution of the FP or CK equations. Similarly, in 
quantum-scale problems, where © (or Q) corresponds to a time-incremental solution for a wave 
function, tp, via a propagator U (see, e.g., [4]) (where ip and U correspond respectively to r/ and p), 
one must typically construct a Hamiltonian appropriate to the quantum system of interest. 

2. Example I: Cole-Hopf/KPZ transformation 

In order to illustrate the utility of ([3D (or @) m solving nonlinear problems, we choose two relatively 
simple embodiments of ([1]), the noise-free Burger's and KPZ equations, and show how the solution, 
([3]), encompasses, and can be used to derive, the CH/KPZ transformations. Importantly, this 
exercise exposes two significant features. First, it illustrates a well known argument [8, 16] for 
obtaining a typical Gaussian transition density via solution of the FP equation; again, determination 
of p, either analytically by solution of the FP or CK equations, or by Monte Carlo simulation 
(leading to comprises the crucial step in solving dU) via ([3]) or ©. 

Second, and reflecting what appears to be a new result, this exercise shows that in some cases, 
one can jump from a time-incremental solution, as in ([3|) or to a non-incremental solution, 
as in the CH/KPZ transform solutions. The idea centers on circumventing a path integral, 
E Ie j+1 u (x{ s ')i s')ds' , over an a priori unknown field, u, by first introducing a transformation 
(j)(x',6) = exp(/(u)) into the representation ([3]), and then by insisting that the stochastic con- 
struction of 4> be purely diffusive, i.e., that the stochastic process, dx^is), used to construct be 
governed by dx<f>{s) = \/2dw(s). 

Thus, considering the noise-free Burger's equation, ue + uu x + u xx = 0, we set rj = u in ([3]), 
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where, for hydro dynamic problems, u is the fluid velocity. The solution in (J3j) then takes the form 



This is a time-incremental, stochastic solution to Burger's equation, and as far as we can determine, 
has not been reported. 

In order to determine p : i) we again assume that the solution point (x, 9) is far removed from 
both 5Q and the final time slice D x {9 = T}, ii) again choose the time increment A9 so that 
the local Reynolds number, or its equivalent, satisfies Ret, = O(l), and iii) set the source term 
/ in dU) to 0. Under these conditions, the transition density, p(x, 9j\x', 6j+i), can be determined 
explicitly using the backward Fokker-Planck equation, pg + u(x',6)p x i + p x >x' = 0, subject to the 
condition p(x, 0j\x', 9j) = S(x — x'), with the result [8,16]: p(x,9j\x' ,9j + \) = [\* / 4irA9}~ 1 exp —[(x' — 
x) 2 / (4A0)], where again, x = x + E Jg j+1 s')ds' is the average position of the RW swarm at 

9j + \. Again, since E fg j+1 u(x(s'), s')ds' = u(x,9j)A9 + 0(A9 2 ), the presence of the path integral 
does not prevent time-incremental solutions; various strategies, involving, e.g., Taylor expansions 
relating the unknown u(x,9j) to known u(x,9j+\), can be easily envisioned. 

Next, arriving at the crucial step in the derivation of the Cole-Hopf solution, we take advantage 
of the general nature of ([3|), and as mentioned above, introduce a transformation 



into ([3]). Again, in order to circumvent the path integral over unknown u, we force the stochastic 




(5) 




(6) 
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construction of <p to be purely diffusive, cpg + <p xx = 0. Importantly, defining a purely diffusive <f> 
of the form given in ([6]) allows us, on one hand, to stretch the solution interval to any arbitrary 
length (by eliminating the path integral), while on the other, providing us with sufficient freedom 
to identify transformations, f(u), between the simple linear heat equation governing eft, and more 
complicated, nonlinear equations governing u. It appears likely that a similar strategy can be used 
to find long-time (non-incremental) transformations and solutions for other nonlinear equations. 
Thus, recognizing that p is now governed by pg + p xx = 0, and inserting <\> in we obtain 



<P(x,6 j ) = I <t>(x',0 j+1 ) yM0j+i-8j)] ^xp-iix' -x) 2 /(A(9 3+1 - 3 ))]dx' (7) 



where now x = x. 

At this stage, several equivalent paths can be taken to complete the derivation. The most 
transparent, which also provides the KPZ transform, first determines f(u) by finding an equivalent 
transformation, eft = <p(f(v)), that solves the related forward-time, noise-free KPZ equation, vt + 
— v xx = 0. This is followed by use of the well-known transformation [17], u = v x , between 
the KPZ and Burger's equations, after which <fi is introduced into ([3]), 9j+i stretched to the final 
backward time T, and the final result obtained. 

Thus, the forward-time form of the governing equation for (j), fa — <j> xx = 0, yields vt — [f v + 
fw/fv]v x — Vxx = 0, which, in turn, transforms to the KPZ equation by choosing [f v +fw/fv] = — 1/2 
(where we again set the coefficient multiplying v xx to 1). Note, we change to a forward time 
solution for <j> in order to obtain the correct signs on v. This is acceptable (and an advantage of 
this approach) since the purely diffusive stochastic process, x<j>i used to construct <p can be run 
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in either the forward or backward time direction. In order to obtain a typical form of the KPZ 
transformation and solution [17], we choose f v = —1/2, which then yields <ft = exp(— v/2), or 
v(x,t) = — 21n</>, identical to the form given in [17] (with v = 1). 

Next, restate eft in terms of u, 4> = ex P [ — 1/2 J x udx~j , and insert this into @. Then iso- 
late u(x,9j) using u(x,9j) = —2(j) x /(f)\ Xj 0., set 6j = t, 6j + i = T, and finally obtain, after some 
manipulation, the classic Cole-Hopf solution of Burger's equation (see, e.g., [12] ): 

fZ^eM-G/2)dx> 

u{x,t) = = (8) 

J_ oo exp(-G/2)(ix / 

where G = G(x'; ,x,t) = u(y, t = 0)dy + ^ 2 t ^ • 
3. Example II: Feynman-Kac formula 

Turning next to solution of nonhomogeneous versions of (H|) via ([3]), we illustrate by using ([3]) 
to derive the Feynman-Kac formula. In particular, consider a one-dimensional Schrodinger-like 
transport equation of the form t]q + bij x + vr] xx = V(x, 6)rj, where b and v either correspond to, or 
are local analogs to drift and diffusivity, and where V is a deterministic function. The Kolmogorov 
solution ^ over A# is thus stated as 



V (x,9 J )= I p(x, 9j\x', 9 j+1 ))r]{x' , 9 j+1 )dx' - E / 3+ V( X (s'), s')r,( X (s% s')ds'. (9) 

-oo 



For present purposes, there is no need to specify p. Thus, we first note that E § e 3+1 V(x(s'), s')r](x(s'), s')ds' 
V(x,9j)r)(x,9j)A9 + 0(A9 2 ). Hence, to 0(A9 2 ), © can be expressed as 

7]{x,9j) = (1 + VA9)- l J^ oo p{x,9 j \x',9 j+l ))i 1 (x\9 j+l )dx', or, using (1 + VA9)- 1 = (1 - VA9) + 
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0{A9 2 ) = exp[-V(x, 9 j )A9]+0{A9 2 ), as n(x, 9j) = exp[-V(x, 9j)A9] fH^pfa 9j\x', 9 j+1 ))r](x', 9 j+1 )dx'+ 
O{A0 2 ). 

Next, taking the expectation and suppressing the error estimate, we rewrite the preceding ex- 
pression for rj(x,0j) as n(x,9j) = exp[— V(x, 0j)A9]r](xj+i, Oj+i), where the average position of the 
RW swarm at Oj+i is expressed as With this last expression, we can now connect rj(x,9j) 

to t](xn,9n), the value of n in the final time-slice D x {9 = T}, evaluated at the swarm's final 
mean position, xn', thus, rj(x,9j) = exp — A9[V(x,9j) + V(xj + i, 9j + \) + V(xj + 2, + • • • + 
V(x N ,9 N )] rj(x N ,6 N ), or finally, rj(x,9j) = E exp [ - J e N V(x{s'), s')ds'] r](x N ,9 N ). 

A closing note regarding nonhomogeneous problems: as before, in cases where an analytical so- 
lution for the transition density, p, cannot be obtained, then an incremental stochastic construction 
can be carried out using an appropriate SDE combined with the approximate representation in ((3J). 
In the present illustration, for example, the stochastic processes used to construct r/ are governed 
by d X {9) = b{ x {9), 9)d9 + s/2vdw(6). 
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